---
title: "Yu, Whang, and Lee_logfile"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
### When Does Audience Matter? Challengers' Stability and Audience Costs
### Replication Main Result

try(setwd(dirname(rstudioapi::getActiveDocumentContext()$path)))
rm(list=ls())

### Load Packages
library(Hmisc)
library(lmtest)
library(sandwich)
library(multiwayvcov)
library(stargazer)
library(jtools)
library(dplyr)
library(tidyr)
library(colorspace)

### Load Dataset
data = read.csv('Yu, Whang, and Lee_dataset.csv')

### Table 1: Correlation between Regime and Stability
rcorr(cbind(data$polity2.1,data$stability.1,data$sfirev.1))

### Table 2: Effect of Stability on Reciprocation

## Without Control
mb1 = glm(recip~poly(stability.1,2)-1
          +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
          +styear_2000+styear_2001+styear_2002+styear_2003
          +styear_2004+styear_2005+styear_2006+styear_2007
          +styear_2008+styear_2009+styear_2010
          ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
mb2 = glm(recip~poly(sfirev.1,2)-1
          +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
          +styear_2000+styear_2001+styear_2002+styear_2003
          +styear_2004+styear_2005+styear_2006+styear_2007
          +styear_2008+styear_2009+styear_2010
          ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

# Robust SEs
mb1.c = coeftest(mb1,vcov=cluster.vcov(mb1, ~ ccode1+ccode2,use_white = T))
mb2.c = coeftest(mb2,vcov=cluster.vcov(mb2, ~ ccode1+ccode2,use_white = T))

## With Control
m1 = glm(recip~poly(stability.1,2)+polity2.1+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m2 = glm(recip~poly(sfirev.1,2)+polity2.1+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

# Robust SEs
m1.c = coeftest(m1,vcov=cluster.vcov(m1, ~ ccode1+ccode2,use_white = T))
m2.c = coeftest(m2,vcov=cluster.vcov(m2, ~ ccode1+ccode2,use_white = T))

## Full Results
stargazer(mb1,m1,mb2,m2,no.space = T
          ,se=list(mb1.c[,2],m1.c[,2],mb2.c[,2],m2.c[,2])
          ,type='text',
          omit='styear',
          order=c(1,2,14,15,3:13))

### Figure 1: Reciprocation ~ Stability
frame1 = data[rownames(model.frame(m1)),]
frame2 = data[rownames(model.frame(m2)),]

pre = make_predictions(m1,pred='stability.1',data=frame1,cluster=c('ccode1','ccode2'))

hist(frame1[frame1$democracy==1,]$stability.1,breaks=20,col='grey',border='white',
     xlim=range(frame1$stability.1),ylim=c(0,50),
     xlab='',ylab='',xaxt='n',yaxt='n',main='')
hist(frame1[frame1$democracy==0,]$stability.1,breaks=20,add=T,col=rgb(1,1,1,alpha=0))
par(new=T)
plot(recip~stability.1,data=pre,type='l',
     ylab='Probability of Target Reciprocation',xlab='Stability',
     xlim=range(frame1$stability.1),ylim=c(0.5,1),lwd=2)
lines(ymax~stability.1,data=pre,lty='dashed',lwd=2)
lines(ymin~stability.1,data=pre,lty='dashed',lwd=2)

legend('topright',legend=c('Democracy','Non-democracy'),fill=c('grey','white'))

### Figure 2: Reciprocation ~ SFI
pre2 = make_predictions(m2,pred='sfirev.1',data=frame2,cluster=c('ccode1','ccode2'))

hist(frame2[frame2$democracy==1,]$sfirev.1,breaks=20,border='white',col='grey',
     xlim=range(frame2$sfirev.1),ylim=c(0,100),
     xlab='',ylab='',xaxt='n',yaxt='n',main='')
hist(frame2[frame2$democracy==0,]$sfirev.1,breaks=20,add=T,col=rgb(1,1,1,alpha=0))
par(new=T)
plot(recip~sfirev.1,data=pre2,type='l',
     ylab='Probability of Target Reciprocation',xlab='State Fragility Index',
     ylim=c(0.5,1),xlim=range(frame2$sfirev.1),lwd=2)
lines(ymax~sfirev.1,data=pre2,lty='dashed',lwd=2)
lines(ymin~sfirev.1,data=pre2,lty='dashed',lwd=2)

legend('topright',legend=c('Democracy','Non-democracy'),fill=c('grey','white'))

### Descriptive Statistics (Full)
data$gdppc = log(data$gdppc)

m.s = glm(recip~sfirev.1+stability.1+polity2.1+
            cinc.ratio+affinity+dist.log+alliance+
            sqevaluation.1+sqevaluation.2+territory+gdpgrowth+
            gdppc+history,data=data,family=binomial(link='probit'))
stargazer(data[!is.na(data$sfirev.1)|!is.na(data$stability.1),
               c('recip',names(m.s$coefficients[-1]))],type='text')

### Table A1: Stability model data
stabilitydata = data[!is.na(data$stability.1),c('recip',names(m.s$coefficients[-1]))]
stargazer(stabilitydata,omit.summary.stat=c('p25','p75'),omit='sfi'
          ,type='text')

length(unique(data[!is.na(data$stability.1),]$ccode1))
length(unique(data[!is.na(data$stability.1),]$ccode2))
range(data[!is.na(data$stability.1),]$styear)

# Table A2: SFI model data
sfidata = data[!is.na(data$sfirev.1),c('recip',names(m.s$coefficients[-1]))]
stargazer(sfidata,omit.summary.stat=c('p25','p75'),omit='stability'
          ,type='text')

length(unique(data[!is.na(data$sfirev.1),]$ccode1))
length(unique(data[!is.na(data$sfirev.1),]$ccode2))
range(data[!is.na(data$sfirev.1),]$styear)
```

```{r}
### Replication Comparison between Linear and Curvilinear Specifications

rm(list=ls())

### Reload Dataset
data = read.csv('Yu, Whang, and Lee_dataset.csv')

### Table OA1: Comparison between Linear and Curvilinear Specifications

## Linear Specification
m1 = glm(recip~stability.1+polity2.1+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m2 = glm(recip~sfirev.1+polity2.1+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m1.c = coeftest(m1,vcov=cluster.vcov(m1, ~ ccode1+ccode2,use_white = T))
m2.c = coeftest(m2,vcov=cluster.vcov(m2, ~ ccode1+ccode2,use_white = T))

## Curvilinear Specification (Main Result)
m3 = glm(recip~poly(stability.1,2)+polity2.1+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m4 = glm(recip~poly(sfirev.1,2)+polity2.1+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m3.c = coeftest(m3,vcov=cluster.vcov(m3, ~ ccode1+ccode2,use_white = T))
m4.c = coeftest(m4,vcov=cluster.vcov(m4, ~ ccode1+ccode2,use_white = T))

## Results
stargazer(m1,m3,m2,m4,no.space=T,
          se=list(m1.c[,2],m3.c[,2],m2.c[,2],m4.c[,2]),
          omit='styear',
          type='text')

```

```{r}
### Replication Stability Trend

rm(list=ls())

### Reload Dataset
data = read.csv('Yu, Whang, and Lee_dataset.csv')

## Threat Variable
threat = data %>%
  select(ccode1,styear) %>%
  distinct() %>%
  mutate(threat = 1)

## Reciprocation Variable
recip = data %>%
  filter(recip == 1) %>%
  select(ccode1,endyear) %>%
  distinct() %>%
  mutate(recip = 1)

## Stability Variable (need to install vdemdata package)
stab = vdemdata::vdem %>%
  select(COWcode,year,e_wbgi_pve) %>%
  drop_na() %>%
  distinct()

## Dataset including Observations without Conflict
trend = data.frame() %>%
  expand(ccode = unique(stab$COWcode),
         year = max(stab$year):min(stab$year)) %>%
  left_join(stab,by=c('ccode'='COWcode','year')) %>%
  left_join(threat,by=c('ccode'='ccode1','year'='styear')) %>%
  left_join(recip,by=c('ccode'='ccode1','year'='endyear')) %>%
  mutate(threat = ifelse(is.na(threat),0,threat)) %>%
  mutate(recip = ifelse(is.na(recip),0,recip)) %>%
  mutate(s = e_wbgi_pve) %>%
  arrange(ccode,year) %>%
  group_by(ccode) %>%
  ## Stability from t-5 to t+5
  mutate(sb1 = lag(e_wbgi_pve,1)) %>%
  mutate(sb2 = lag(e_wbgi_pve,2)) %>%
  mutate(sb3 = lag(e_wbgi_pve,3)) %>%
  mutate(sb4 = lag(e_wbgi_pve,4)) %>%
  mutate(sb5 = lag(e_wbgi_pve,5)) %>%
  mutate(sa1 = lead(e_wbgi_pve,1)) %>%
  mutate(sa2 = lead(e_wbgi_pve,2)) %>%
  mutate(sa3 = lead(e_wbgi_pve,3)) %>%
  mutate(sa4 = lead(e_wbgi_pve,4)) %>%
  mutate(sa5 = lead(e_wbgi_pve,5)) %>%
  ungroup()

s = seq(min(trend$s,na.rm=T),max(trend$s,na.rm=T),len=11)
t = 1:11

## Mean of Stability before/after Threat (Reciprocation)
threat.agg = apply(apply(trend[trend$threat==1,],1,
                         function(x) x[c(11:6,12:16)]),1,mean,na.rm=T)
recip.agg = apply(apply(trend[trend$recip==1,],1,
                        function(x) x[c(11:6,12:16)]),1,mean,na.rm=T)

### Figure OA1

### Panel A: Variation of Stability Before/After Threat (N=75)

plot(s~t,type='n',
     ylab='Stability',
     xlab='',xaxt='n',
     main='A. Before/After Threat (N=75)')
for(i in sample(1:nrow(trend[trend$threat==1,]),75,replace=F)){
  ss = trend[i,c(paste('sb',5:1,sep=''),
                 's',
                 paste('sa',1:5,sep=''))]
  lines(as.vector(as.matrix(ss))~t,
        col=adjust_transparency('gray',alpha=0.5))  
}
abline(v=6,col='red',lty='dashed')
lines(threat.agg)
axis(side=1,
     at=1:11,
     labels=c('t-5',
              't-4',
              't-3',
              't-2',
              't-1',
              't',
              't+1',
              't+2',
              't+3',
              't+4',
              't+5'))

### Panel B: Variation of Stability Before/After Reciprocation (N=75)

plot(s~t,type='n',
     ylab='Stability',
     xlab='',xaxt='n',
     main='B. Before/After Reciprocation (N=75)')
for(i in sample(1:nrow(trend[trend$recip==1,]),75,replace=F)){
  ss = trend[i,c(paste('sb',5:1,sep=''),
                 's',
                 paste('sa',1:5,sep=''))]
  lines(as.vector(as.matrix(ss))~t,
        col=adjust_transparency('gray',alpha=0.5))  
}
abline(v=6,col='red',lty='dashed')
lines(recip.agg)
axis(side=1,
     at=1:11,
     labels=c('t-5',
              't-4',
              't-3',
              't-2',
              't-1',
              't',
              't+1',
              't+2',
              't+3',
              't+4',
              't+5'))

### Panel C: Variation of Stability Before/After Threat (N=150)

plot(s~t,type='n',
     ylab='Stability',
     xlab='',xaxt='n',
     main='C. Before/After Threat (N=150)')
for(i in sample(1:nrow(trend[trend$threat==1,]),150,replace=F)){
  ss = trend[i,c(paste('sb',5:1,sep=''),
                 's',
                 paste('sa',1:5,sep=''))]
  lines(as.vector(as.matrix(ss))~t,
        col=adjust_transparency('gray',alpha=0.5))  
}
abline(v=6,col='red',lty='dashed')
lines(threat.agg)
axis(side=1,
     at=1:11,
     labels=c('t-5',
              't-4',
              't-3',
              't-2',
              't-1',
              't',
              't+1',
              't+2',
              't+3',
              't+4',
              't+5'))

### Panel D: Variation of Stability Before/After Reciprocation (N=150)

plot(s~t,type='n',
     ylab='Stability',
     xlab='',xaxt='n',
     main='D. Before/After Reciprocation (N=150)')
for(i in sample(1:nrow(trend[trend$recip==1,]),150,replace=F)){
  ss = trend[i,c(paste('sb',5:1,sep=''),
                 's',
                 paste('sa',1:5,sep=''))]
  lines(as.vector(as.matrix(ss))~t,
        col=adjust_transparency('gray',alpha=0.5))  
}
abline(v=6,col='red',lty='dashed')
lines(recip.agg)
axis(side=1,
     at=1:11,
     labels=c('t-5',
              't-4',
              't-3',
              't-2',
              't-1',
              't',
              't+1',
              't+2',
              't+3',
              't+4',
              't+5'))

```

```{r}
### Replication Democracy Variables

rm(list=ls())

### Reload Dataset
data = read.csv('Yu, Whang, and Lee_dataset.csv')

### Polity IV source

## democracy (Polity2 >= 7)

m11 = glm(recip~poly(stability.1,2)+democracy+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m12 = glm(recip~poly(sfirev.1,2)+democracy+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m11.c = coeftest(m11,vcov=cluster.vcov(m11, ~ ccode1+ccode2,use_white = T))
m12.c = coeftest(m12,vcov=cluster.vcov(m12, ~ ccode1+ccode2,use_white = T))

## exrec (Executive Recruitment)

m13 = glm(recip~poly(stability.1,2)+exrec+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m14 = glm(recip~poly(sfirev.1,2)+exrec+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m13.c = coeftest(m13,vcov=cluster.vcov(m13, ~ ccode1+ccode2,use_white = T))
m14.c = coeftest(m14,vcov=cluster.vcov(m14, ~ ccode1+ccode2,use_white = T))

## exconst (Executive Constraints)

m15 = glm(recip~poly(stability.1,2)+exconst+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m16 = glm(recip~poly(sfirev.1,2)+exconst+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m15.c = coeftest(m15,vcov=cluster.vcov(m15, ~ ccode1+ccode2,use_white = T))
m16.c = coeftest(m16,vcov=cluster.vcov(m16, ~ ccode1+ccode2,use_white = T))

## polcomp (Political Competition and Opposition)

m17 = glm(recip~poly(stability.1,2)+polcomp+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m18 = glm(recip~poly(sfirev.1,2)+polcomp+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m17.c = coeftest(m17,vcov=cluster.vcov(m17, ~ ccode1+ccode2,use_white = T))
m18.c = coeftest(m18,vcov=cluster.vcov(m18, ~ ccode1+ccode2,use_white = T))

### V-Dem source

## v2x_polyarchy (Electoral)

m21 = glm(recip~poly(stability.1,2)+v2x_polyarchy+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m22 = glm(recip~poly(sfirev.1,2)+v2x_polyarchy+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m21.c = coeftest(m21,vcov=cluster.vcov(m21, ~ ccode1+ccode2,use_white = T))
m22.c = coeftest(m22,vcov=cluster.vcov(m22, ~ ccode1+ccode2,use_white = T))

## v2x_libdem (Liberal)

m23 = glm(recip~poly(stability.1,2)+v2x_libdem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m24 = glm(recip~poly(sfirev.1,2)+v2x_libdem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m23.c = coeftest(m23,vcov=cluster.vcov(m23, ~ ccode1+ccode2,use_white = T))
m24.c = coeftest(m24,vcov=cluster.vcov(m24, ~ ccode1+ccode2,use_white = T))

## v2x_partipdem (Participatory)

m25 = glm(recip~poly(stability.1,2)+v2x_partipdem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m26 = glm(recip~poly(sfirev.1,2)+v2x_partipdem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m25.c = coeftest(m25,vcov=cluster.vcov(m25, ~ ccode1+ccode2,use_white = T))
m26.c = coeftest(m26,vcov=cluster.vcov(m26, ~ ccode1+ccode2,use_white = T))

## v2x_delibdem (Deliberate)

m27 = glm(recip~poly(stability.1,2)+v2x_delibdem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m28 = glm(recip~poly(sfirev.1,2)+v2x_delibdem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m27.c = coeftest(m27,vcov=cluster.vcov(m27, ~ ccode1+ccode2,use_white = T))
m28.c = coeftest(m28,vcov=cluster.vcov(m28, ~ ccode1+ccode2,use_white = T))

## v2x_egaldem (Egalitarian)

m29 = glm(recip~poly(stability.1,2)+v2x_egaldem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m30 = glm(recip~poly(sfirev.1,2)+v2x_egaldem+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m29.c = coeftest(m29,vcov=cluster.vcov(m29, ~ ccode1+ccode2,use_white = T))
m30.c = coeftest(m30,vcov=cluster.vcov(m30, ~ ccode1+ccode2,use_white = T))

### Boix et al. (e_boix_regime)

m31 = glm(recip~poly(stability.1,2)+e_boix_regime+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$stability.1),],family=binomial(link='probit'))
m32 = glm(recip~poly(sfirev.1,2)+e_boix_regime+cinc.ratio+affinity+dist.log+alliance
         +sqevaluation.1+sqevaluation.2+territory+gdpgrowth
         +log(gdppc)+history-1
         +styear_1995+styear_1996+styear_1997+styear_1998+styear_1999
         +styear_2000+styear_2001+styear_2002+styear_2003
         +styear_2004+styear_2005+styear_2006+styear_2007
         +styear_2008+styear_2009+styear_2010
         ,data=data[!is.na(data$sfirev.1),],family=binomial(link='probit'))

m31.c = coeftest(m31,vcov=cluster.vcov(m31, ~ ccode1+ccode2,use_white = T))
m32.c = coeftest(m32,vcov=cluster.vcov(m32, ~ ccode1+ccode2,use_white = T))

### Table OA4
stargazer(m11,m13,m15,m17,m31,omit='styear',
          se=list(m11.c[,2],
                  m13.c[,2],
                  m15.c[,2],
                  m17.c[,2],
                  m31.c[,2]),
          type='text',
          no.space = T)

### Table OA5
stargazer(m12,m14,m16,m18,m32,omit='styear',
          se=list(m12.c[,2],
                  m14.c[,2],
                  m16.c[,2],
                  m18.c[,2],
                  m32.c[,2]),
          type='text',
          no.space = T)

### Table OA6
stargazer(m21,m23,m25,m27,m29,omit='styear',
          se=list(m21.c[,2],
                  m23.c[,2],
                  m25.c[,2],
                  m27.c[,2],
                  m29.c[,2]),
          type='text',
          no.space = T)

### Table OA7
stargazer(m22,m24,m26,m28,m30,omit='styear',
          se=list(m22.c[,2],
                  m24.c[,2],
                  m26.c[,2],
                  m28.c[,2],
                  m30.c[,2]),
          type='text',
          no.space = T)

## Descriptive Statistics (Full)
data$gdppc = log(data$gdppc)

m.s = glm(recip~sfirev.1+stability.1+
            democracy+exrec+exconst+polcomp+e_boix_regime+
            v2x_polyarchy+v2x_libdem+v2x_partipdem+v2x_delibdem+v2x_egaldem
            +cinc.ratio+affinity+dist.log+alliance+
            sqevaluation.1+sqevaluation.2+territory+gdpgrowth+gdppc+history,
          data=data,family=binomial(link='probit'))
stargazer(data[!is.na(data$sfirev.1)|!is.na(data$stability.1),
               c('recip',names(m.s$coefficients[-1]))],type='text')

# Table OA2: Stability model data
stabilitydata = data[!is.na(data$stability.1),c('recip',names(m.s$coefficients[-1]))]
stargazer(stabilitydata,omit.summary.stat=c('p25','p75'),omit='sfi'
          ,type='text',no.space = T)

length(unique(data[!is.na(data$stability.1),]$ccode1))
length(unique(data[!is.na(data$stability.1),]$ccode2))
range(data[!is.na(data$stability.1),]$styear)

# Table OA3: SFI model data
sfidata = data[!is.na(data$sfirev.1),c('recip',names(m.s$coefficients[-1]))]
stargazer(sfidata,omit.summary.stat=c('p25','p75'),omit='stability'
          ,type='text',no.space = T)

length(unique(data[!is.na(data$sfirev.1),]$ccode1))
length(unique(data[!is.na(data$sfirev.1),]$ccode2))
range(data[!is.na(data$sfirev.1),]$styear)

```